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C3 . Abstract 



We present a mathematical model for solvated crystallization of a-lactose monohydrate 
in semi-batch mode. The process dynamics are governed by conservation laws including 
population, molar and energy balance equations. We present and discuss the model and then 
control the process with the goal to privilege the production of small particles in the range 
between 10 -5 and W~ 4 fim. We compare several specific and unspecific cost functions leading 
to optimal strategies with significantly different effects on product quality. Control inputs 
are temperature, feed rate, and the choice of an appropriate crystal seed. 
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1 Introduction 



in 

o: 

CO ■ Crystallization is the unitary operation of formation of solid crystals from a solution. In process 
\ engineering crystallization is an important separation process used in chemical, pharmaceutical, 

•l-H ■ 

^ ■ food, material and semiconductor industries. Mathematical models are described by conserva- 

h: 

5r ■ tion laws with population, molar and energy balance equations. Crystallizers can be operated 
either in batch, semi-batch or continuous mode. Semi-batch crystallization is widely used in the 
pharmaceutical and fine chemical industry for the production of solids in a variety of operating 
modes. 

Crystallization processes are described by balance equations, including a population balance 
for the particle size distribution, a molar balance for the distribution of solute, and an energy 
balance equation to model thermodynamic phenomena. In the food-processing industry, there has 
been a growing interest in the crystallization of lactose in recent years [4,5,7]. In this paper we 
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study a model of solvated crystallization of a-lactose monohydrate, which includes four interacting 
populations, one of them aging, in tandem with an energy balance. Two forms of lactose (a- and 
/3-lactose) exist simultaneously in aqueous solution, the exchange being governed by mutarotation. 
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Figure 1: Solvated crystallization of a-lactose monohydrate with complex population dynamics 
featuring one aging and three ageless populations. Exchange rates depend on temperature. Con- 
trols are feed rate in semi-batch mode, temperature of envelope, and the distribution of the crystal 
seed. Based on mathematical modeling the process is optimized in order to maximize the particle 
mass in a given small size range L\ ow < L < Lhigh- 



Crystallization and precipitation processes are modeled as highly nonlinear and complex dy- 
namical systems. This makes it interesting to simulate, control and optimize these processes in 
order to enhance product quality in various situations [1,2]. In order to control crystallization 
processes, Hu and Rohani [3] have studied different heuristic cooling methods such as linear cool- 
ing, natural cooling and controlled cooling. In this study we control of the process of formation 
of a-lactose crystals in such a way that the growth of very large crystals is avoided and the bulk 
of crystal mass occurs in a small particle range between and 10 -4 /im. In order to achieve 
this goal we use optimal control of the process in semi-batch mode by variations of temperature, 
feed rate, and also by an appropriate choice of the crystal seed, and based on a variety of different 
figure of merit functions. 



2 Dynamic model of process 

In this section the dynamic model of semi-batch crystallization of a-lactose monohydrate is de- 
scribed by presenting the population, molar and energy balance equations. 
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Population balance equation: The population balance equation is a first-order PDE 

| (V(t)n(L, t)) + V(t)G (c a (t), cp{t),T(t)) = V(t)n(L, f)± (1) 

where n(L,t) is the distribution of a-lactose crystals, c a (t),cp(t) are the dimensionless concentra- 
tions of a- and /3-lactose in the liquid phase, V(t) is the volume of slurry in the crystallizer, a 
dependent variable given in (11), G (c a ,cp,T) is the temperature-dependent growth coefficient of 
a-crystals, assumed independent of crystal size L, and the right hand side represents source and 
sink terms. We add the boundary condition 

_ B(c a (t),cp(t),T(t)) 

n ^ t} - G(c a (t),cp(t),T(t)Y {I) 

and the initial condition 

n(L,0)=n (L), (3) 

where n (L) is the crystal seed. It is convenient to introduce the moments of the crystal size 
distribution function 

POO 

fi v (t)= / n(L,t)L u dL, i/ = 0,l,..., 
Jo 

which allows to break (1) into an infinite sequence of ODEs if the source and sink terms may be 
neglected. In the present situation this amounts to neglecting agglomeration and also attrition 
effects. Using (2) we obtain the equations 

d ^ + V§)^ {t) " uG (Ca(t) ' C/?(t) ' T{t)) ^- l(t) = °' (4) 

v — 1, 2, . . . , in tandem with 

dtio(t) , V\t) 



dt V(t) 
The initial conditions are then 



fi o {t)-B{c a {t),cp{t),T{t)) = 0. (5) 



POO 

H v (0)= / n (L)L u dL, v = 0, 1, . . . . 
Jo 



Solvent mass balance: 

dm B2 o(t) 



{R- 1 - l)3k vPcry G (c a (t), c p (t),T(t)) V(t)n 2 (t) + q R2 o{t) (6) 



dt 

Here mH 2 o is the mass of water in the aqueous solution, which changes due to feed, qn 2 o, and due 
to the integration of water molecules into the a-crystals. The constant R = M CTy /M a = 1.0525 is 
the ratio of the molar masses of the solid and liquid phases of a-lactose. 



Concentration of a-lactose: The dimensionless concentration of a-lactose c a in the solution is 



defined as m a = c a mu 2 o and satisfies the differential equation 

= -^-fcMW ~ R' 1 ) ~ R- 1 }^^ - h(T(t))c a (t) + k 2 {T{t))c,{t) (7) 
at rnu 2 o\t) at 

+ (c iW - c «w) 9Hao(t) 



Here m cry is the crystal mass in the slurry, c+ is the feed rate of a lactose during the semi-batch 
phase, and ki(T), k 2 (T) are the temperature dependent mutarotation exchange rates between a- 
and /3-lactose in the liquid phase. 

Concentration of /3-lactose: The dimensionless concentration of /3-lactose eg is defined as mp = 
c /3 m n 2 o an d satisfies the differential equation 

«h 2 o(*) 



+ (c+(o - c „(f)) 



m Ha o(*) 

Energy balance: The temperature hold system describes the interaction between crystallizer 
temperature, the temperature of the jacket, and the control signal, to which the internal heat 
balance due to enthalpy is added. We have 



dT{t) = m 



P 2 (t)(T(t) - T rcf ) - AH^^- + UA(t) (T jackct (t) - T(t)) (9) 



dt 

+ te 2 o(i) (Cg 2 o + C p a c a (0) + CJc^(O)) (T fccd - T ref ) 
where 

rfTja g t(t) = -0.0019(T jackct (t) - T sp (t)) (10) 

was obtained through identification of the system. Here T(t) is the temperature of the slurry, T ref = 
25°C a constant reference temperature, Tf ee d the temperature of feed, which is the temperature of 
H 2 in this case, assumed constant in this study, Tj ackct (t) is the temperature of the crystallizer 
jacket, and T sp (t) is the set point temperature, which is used as a control input to regulate Tj acket (t), 
and therefore indirectly T(t) via the heat exchange between the envelope and the crystallizer 
through the contact surface A(t), which is determined through V(t). The constants C^ 2Q , C£, Cp 
are the specific heat capacities. We have 

P^t)- 1 = m H2 o(t)C£ 20 + m a (t)C p a + mp(t)C p + m CIy (t)C p cvy , 
_ dm H2 o(t) dm a (t) dm^jt) dm cry {t) 
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with m a = c Q m H2 o, rrip = c^m^o 

Mutarotation: The mutarotation exchange coefficients hi, h 2 are temperature dependent and are 
determined as 

h 2 (T) = h ■ exp(-E/(R ■ (T + 273.15))), 

h m (T) = 1.64 - 0.0027 • T, h(T) = h 2 {T) ■ h m (T). 
The equilibrium of mutarotation therefore occurs at 

10.9109 • exp(0.02804 • T) 
c Q)Sa t, eq (T) - 100(1 + h m {T)) ' 

F(T) = 0.0187 • exp(0.0236 • T), 

c a , sa t(c/3,T) = c aiSat , cq (T) - F(T)(cp - h m (T) 

X Q*,sat ,eq ( T^j ) . 

Nucleation rate: The nucleation rate is based on a phenomenological law 

B 



B(c a ,cp,T) = h b exp 



(T + 273.15)3 In 2 (^^) 



as is the growth rate 
Growth rate: 



G(c a ,cp,T) = h g (c a - c ajSa ,t(cp,T)) . 



Volume: The total volume of slurry V(t) is a dependent variable, which can be expressed as a 
function of the states c a , cp and m H2 o through 

Therefore 

^ = 3k v G(ca(t),cp(t),T(t))V(t)n 2 (t) + [p-^ aCa (t) +p^c p (t) + p£ Q ] 

_! dc a (t) _! dcp(t)\ 



/ \ ( i dco^itj i 



dt J ' 

Crystal mass: The total crystal mass satisfies the equation 

3h vPcTy G(c a (t), cp(t),T(t))V(t)^ 2 {t). (12) 



dm CTy (t) 



dt 
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Table 1: Numerical constants 



3 Optimal control problem 

The benefit of the moment approach is that we may choose a finite number of moment equation to 
replace (1). Our present approach is to retain a sufficient number of moments so that the salient 
features of the seed n (L) may be captured by these moments, and in our experiments we decided 



6 



n(L t) 


H / 77? 777^ 


n^rtirlp sit;p Hi < ^t,ri"hnt,inn 


m a (t) 


kg 


mass of a-lactose in solution 




kg 


mass of /3-lactose in solution 


V(t) 


kg 


volume of slurry 


A(t) 


9 

m 


contact surface 



TABLE 2: Units of dynamic quantities 



to retain the moments /x , • • • , A*5- The remaining states of the system dynamics are then m H2 o, 
m CTy , c a , c^, T, T jacket . The control inputs are u x = T sp and u 2 = <?h 2 o 
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Table 3: Initial values for study 1 
In this work, we compare between several policies : 

• Policy 1 : We fix the values of the set-point temperature T sp = 15 [°C] and the feed rate of 
solvent qn 2 o = 0.0056[Kg/h]. This policy is referred to as constant in the figures. 

• Policy 2 : Here we fix the value of the feed rate of solvent gn 2 o = 0.0056 [Kg/h], while the 
set-point temperature T sp (t) starts at T sp (0) = 15 [°C] and decreases linearly. This policy is 
called linear in the figures. 

• Policy 3 : We control the set-point temperature u±(t) = T sp (t) and also the feed rate of 
solvent u 2 (t) = <?h 2 o(£) using various objectives. This policy is called optimal in the figures. 
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3.1 Scenario 1 

Our first control problem minimizes the weighted mean size diameter c/43 = — at fixed final time 
tf — 11000 seconds. This is cast as the optimization program 

/j 4 (t/) 



minimize 



d A3 (t f ) 



subject to dynamics (4) — (12) 
< V{t) < V max 

0°C < T(t) < 70°C ( 13 ) 

(*)> Comsat 

0°C < T sp (t) < 40°C 
< qn 2 o(t) < 0.1 

The control variables are set-point temperature U\(t) = T sp (t) and feed rate of water U2(t) = 
QK 2 o(t). The percentages c+ and c^" of lactose in the feed are kept constant. 

3.2 Scenario 2 

Our second control problem minimizes the nucleation rate B(c a ,C/3,T) at the fixed final time 
tf = 11000 seconds. This is cast as the optimization program 



minimize Bit A 

V /7 (14) 

subject to constraints of (13) 
The control variables are again T sp and gH 2 o 



3.3 Scenario 3 

Our third control problem minimizes the coefficient of variation CV at the fixed final time tf = 
11000 seconds. This is the optimization program 

minimize CV(t f ) = ^M^M - 1 

1 /; (M*/)) 2 (15) 

subject to constraints of (13) 

The control variables are again T sp and gH 2 o- 

In Figures 2, 3 we present results obtained with the optimal regulation of set-point temperature 
Ui = T sp and feed rate u 2 = qn 2 o an d compare these to more standard scenarios, where temperature 
and feed rate are fixed or follow simple heuristic profiles proposed in the literature. 
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Figure 2: Optimal set-point temperature profile U\ = T sp for minimization of the three criteria 
B, CV and 0^43 with fixed final time tf. 
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Figure 3: Optimal feed profile u<i = qn 2 o f° r minimization of the three criteria B, d^ and CV 



with fixed final time t 



In Figure 4 we present the optimal crystal size distribution obtained from minimization of the 
weighted mean size diameter in (13) compared with standard scenarios. Figure 5 shows the 
crystal size distribution for the optimal control of nucleation rate B and coefficient of variation 
CV compared to the more standard scenarios. The optimal profile for the nucleation rate shows 
the existence of two peaks which indicates the existence of two populations of crystals. 

In Figure 6 we present the evolution of solubility of a— lactose and the temperature of the 
crystallizer by comparing several scenarios. The profile of solubility shows an early peak, which 
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L (n m) L fti m) 

Figure 4: Final crystal size distribution L i— > n(L,tf) displayed for minimization of B (blue) and 
CV (red) compared with linear (magenta) and constant (dashed black) policies for fixed final time 
^finai- Right hand image shows zoom on range [10 2,1 , 10 2 ' 6 ]. 
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Figure 5: Final crystal size distribution L \— > n(L,tf) for minimizing (black solid), compared 
with linear (magenta) and constant (black dashed) policies for fixed final time tfi na i- Right hand 
image shows zoom on range [10 2 * 1 , 10 2 ' 6 ]. 

correspond to a sharp decrease in the temperature profile of the crystallizer. Comparison between 
the cost functions shows that the highest peak occurs when minimizing the weighted mean size 
diameter g? 43 . In the case of minimization of nucleation rate B, we see the existence of two peaks 
which correspond with two peaks on crystal size distribution profile. 

In Figure 7 we present the evolution of nucleation rate B and growth rate G in comparison 
between the several scenarios. At the beginning of the profile of nucleation rate B, we note the 
highest peak in case of minimization of weighted mean size diameter d^. 

In Figure 8 we present the coefficient of variation CV and volume of crystallizer in comparison 
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Figure 6: Left image compares solubility, right image compares temperature of crystallizer for 
the five policies constant (dashed black), linear (magenta), optimal with B (blue), optimal with 
CV (red) and optimal with d i3 (black continuous) for t fina j fixed. 
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Figure 7: Left image compares growth rate, right image compares crystal mass for the five 
policies constant (dashed black), linear (magenta), optimal with B (blue), optimal with CV (red), 
optimal with <i 43 (black continuous) for t final fixed. 

between all objectives and scenarios. The volume profiles show that optimization of different cost 
functions may lead to fairly different ways of filling the crystallizer in the semi-batch phase. For 
instance, filling in the linear policy occurs much faster than for minimization of the coefficient of 
variation, which gives the slowest filling. 

Figure 9 compares the profiles of overall crystals mass and of the weighted mean size diameter 
g?4 3 in all scenarios. 
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Figure 8: Evolution of coefficient of variation CV(t) left, evolution of volume V(t) right. Com- 
parison of the five policies constant (dashed black), linear (magenta), optimal for B (blue), optimal 
for CV (red), optimal for ^43 (black continuous), for t^ na \ fixed. 





Figure 9: Evolution of weighted mean size diameter (i 43 left, evolution of crystal mass right. 
Comparison of the five policies constant (dashed black), linear (magenta), optimal for B (blue), 
optimal for CV (red), optimal for (black continuous), for tfi na j fixed. 



3.4 Scenario 4 

The next extension is to add the moments of Uq{L) as unknown parameters, the idea being that 
a suitable choice of the initial seed of given mass should give even better results. We decide to fix 
the total mass of crystal seed as k v Vop crj J °° n Q (L)L 3 dL = 0.1kg. That leads to the optimization 
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program 

minimize d 43 

subject to constraints of (14) (16) 

POO 

fi 3 (Q) = / n (L)L 3 dL = 0.0812 
Jo 

where now T re f(t), qu 2 o(t) and /io(0), /ii(0), /i2(0), /^(O), /is(0) are optimization variables. 

The interesting point of this program is that once the optimal solution (T* ef , q^ 20 , ij* v ) is reached, 
we need to reconstruct a function Uq(L) such that its moments 0, . . . , 5 coincide with /ig, /4, /ig = 
0.0812, nl, nl. This can be achieved by solving the maximum entropy function reconstruction 
problem 

POO 

minimize / n (L) log no (L) dL 

J °poo (17) 
subject to / L u n (L)dL = /x*, v = 0, . . . , 5. 

Jo 

Notice that (17) may be solved by standard software, see e.g. Borwein et al. [8], [9]. 
3.5 Scenario 5 

The natural figure of merit to maximize the crystal mass within a certain range L\ < L < L2 of 
small particle sizes is 



max / L 3 n(L,tf)dL 

J Li 



at the final time £/, but this objective is not directly accessible in the moment approach. Substrates 
like B, CV, d i3 are non-specific and must be expected to give only a crude approximation of (18). 
We therefore propose the following more sophisticated strategy, which is compatible with the 
moment approach. 

We define a target particle size distribution ni(L), which has a bulk in the range [L 1; L 2 ], 
normalized to satisfy 

POO 

v z = \ L 3 ni (L)dL = 1. 
Jo 

Then we compute as many of its moments uo, . . . ,Un as we wish to use, where as before N = 5 in 
our tests. The optimization program we now solve is 

N 

minimize 2J w i {Vijtf) - /i 3 (t/)z^) 2 

, s (19) 
subject to constraints of (13) 

tf — ^max 

13 



where the Wi are suitably chosen weights. Notice that the least squares objective of (19) tries to 
bring the moments of n(L,tf) as close as possible to the moments of the function fi^(tf)ni(L), 
which has the correct shape, and has the same total crystal mass as n(L,tf). Here the final time 
is considered free. 
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Figure 10: Scenario 5: Different initial seeds Uq(L) shown on the left lead to different products 
n(L,tf) on the right. Blue uses 1x3(0) = 0.291, magenta uses 113(0) = 0.401. 



4 Method 

For our testing we have used the solver ACADO [10] based on a semi-direct single or multiple- 
shooting strategy, including automatic differentiation, based ultimately on the semi-direct multiple- 
shooting algorithm of Bock and Pitt [11]. ACADO is a self-contained public domain software 
environment written in C++ for automatic control and dynamic optimization. 

Alternatively, we also use the solver PSOPT [12], which is a public domain extension of the 
NLP-slover IPOPT [13] or SNOPT [14] and is based on pseudospectral optimization which uses 
Legendre or Cheybyshev polynomials and discretization based on Gauss-Lobatto nodes. 

A difficulty with both solvers is the strong dependence of convergence and solutions on the 
initial guess, as must be expected in a local optimization context. Our testing shows that it is 
often mandatory to have a simulated study (x- m it,u- m i t ) available to start the optimization from 
that point. This initial guess may use parameters from a previous optimization study, which give 
already a decent cost in the present study. In some cases homotopy techniques, using for instance 
tf as a parameter, have to be used. 

Once optimal controls u* = (T* ef , q# 2 o) have been computed in any one of the scenarios, we use 
the full crystallizer model (1), (2), (6) - (10) to simulate the system, using an initial seed n Q (L) 

14 



which produces the initial moments /i„(0). In those cases where the moments of the initial seed 
are parameters, which are also optimized, we use the optimal //* = . . . , fi* N ) to compute an 
estimation Uq{L) of the optimal crystal seed with these moments using [8] and [9]. 

The final stage in each experiment is a simulation of the full population balance model using 
the optimal (T s * p , q^ 2Q , n (L)*), obtained from the moment-based optimal control problem. 

5 Conclusion 

We have presented and tested several control strategies which allow to maximize the crystal mass 
of particles of small size, typically in a range of 10~ 5 — 10~ 4 /im. Our approach was compared to 
more standard heuristic control policies used in the literature to regulate temperature and feed rate 
in semi-batch mode. Our simulated numerical results show that it is beneficial to apply optimal 
control strategies in semi-batch solvated crystallization of a-lactose monohydrate. 
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